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ABSTRACT 

Many flexible parameterizations exist to represent data on the sphere. In addition to the venerable spherical 
harmonics, we have the Slepian basis, harmonic splines, wavelets and wavelet-like Slepian frames. In this paper 
we focus on the latter two: spherical wavelets developed for geophysical applications on the cubed sphere, 
and the Slepian "tree", a new construction that combines a quadratic concentration measure with wavelet-like 
multiresolution. We discuss the basic features of these mathematical tools, and illustrate their applicability in 
parameterizing large-scale global geophysical (inverse) problems. 
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1. INTRODUCTION 

This paper is about parameterization, and its role in geophysical inverse problems: the analysis and representation 
of volumetric properties, with a particular emphasis on the three-dimensional ball and its surface, the two- 
dimensional sphere. This is not the sole purview of geophysics (e.g. geodesy, geodynamics, seismology): there 
is a large literature on the subject in virtually every area of scientific inquiry (e.g. medical imaging, astronomy, 
cosmology, computer graphics, image processing, ...). For this reason we limit ourselves here to a few technical 
aspects, some of which have been published before (refs[l][2| but are illustrated with different examples. 

We begin with a new class of spherical wavelet basis on the ball designed for the analysis of geophysical 
models and for the tomographic inversion of global seismic dataHIHEl Its multiresolution character allows for 
modeling with an effective spatial resolution that varies with position within the Earth. We introduce two types of 
discrete wavelet transforms in the angular dimension of the "cubed sphere" , a well-known Cartesian-to-spherical 
mapping.'^'^ These are applied to analyze the information of terrestrial topography data and in seismic wavespeed 
models of the Earth's mantle across scale space. The localization and sparsity properties of the wavelet bases 
allow finding a sparse solution to inverse problems by iterative minimization of a combination of the £2 norm 
of the data residuals and the £1 norm of the model wavelet coefficients. These have now been validated in 
realistic synthetic experiments to likely yield important gains in the inversion of seismic data for global seismic 
tomography, the procedure by which seismic waveforms are being used to image the three-dimensional distribution 
of elastic (compressional or P and shear or S) wave speeds {Vp and Vs) inside the Earth. 

A second class of transforms is inspired by the type known as Slepian functions. These we understand to 
be families of orthogonal functions that are all defined on a common, e.g. geographical, domain, where they 
are either optimally quadratically concentrated or within which they are exactly limited, and which at the same 
time are exactly confined within a certain bandwidth, or maximally concentrated therein. Originally developed 
for the study of time series by Slepian, Landau and PollalEEl in the 1960s, they have now been fully extended 
to the multidimensional Cartesian plane^ and the surface of the sphere, f^^Jfl^ where most of the geophysical 
applications lie. We have reported on some aspects of spherical Slepian functions and their role in the analysis 
and representation of geophysical data in previous papers in this series EHESl In this contribution we announce 
a hybrid construction that blends desirable aspects and properties of Slepian functions with ideas from wavelet 
and frame theory. We summarize the main ideas contained in Chapter 7, Multiscale Dictionaries of Slepian 
Functions on the Sphere^ of the Ph. D. thesis by Eugene Brevdo,!^ deferring a full publication to a later date. 



2. WAVELETS ON THE CUBED SPHERE 



The use of wavelets is still no matter of routine in global geophysics, beyond applications in one and two Cartesian 
dimensions. This despite there being a wealth of available constructions on the sphereJ^^"^ However, the cited 
studies are mostly concerned with surfaces, not volumes, and many of them require custom design and special 
algorithms. For the application to global seismic tomography that we envisage, we take advantage of the ease and 
flexibility of existing Cartesian constructions by implementing standard algorithms on a spherical-to- Cartesian 
map proposed by Ronchi et al. (1996), a grid that has long since proven its utility in the geosciences and beyond. 

As in ref. |5j we define a coordinate 4-tuple (^, 77, r, for each of the = 1 ^ 6 "chunks" comprising this 
"cubed sphere" at radius r. In ref. 1 we list explicit formulas for the forward and inverse mapping of cubed-sphere 
to Cartesian coordinates. Here we suffice to say that a "master" surface chunk defined by the 2x2^ linearly 
spaced coordinate pairs — 7r/4 < (^,/?) < 7r/4 maps, in a general sense, to the Cartesian coordinate vector 

[x,y,z] = [l tan 77 - tan^] (1 + tan^ ^ + tan^ t^)"^/^ (1) 

which is rotated to occupy a total of six non-overlapping patches tiling the sphere, after which the entire con- 
struction is rotated over a standard set of Euler angles a = 0.0339, /3 = 1.1705, and 7 = 1.1909. This results in 
the configuration shown in Figure [l] (/e/i^). 

Alternatively, for reasons that will be made clear, we start from the coordinates —Stt/S < (^,77) < 37r/8 and 
apply the same rotations, thereby producing six overlapping "superchunks" , as in Figure [l] (n^/i^). Throughout 
this paper we will quote TV as the angular resolution level of our cubed sphere, which implies that it has 6 x 2^^ 
such elements, with typical seismic tomography grids having N = 7. The Euler angles used in our construction 
were chosen for geographical convenience, as can be seen from the continental outlines in Figure |2j 

On these two grids we apply standard Cartesian transforms, e.g. orthogonal^^ and bi-orthogonal^^ ones, 
whereby for the case of the non-overlapping "chunks" construction we accommodate the presence of the "seams' 
by switching to special boundary filters at each of the edges, and applying preconditioners to the data prior 
to transformation in order to guarantee the usual polynomial cancellation throughout the closed rectangular 
interval, as is well established.^ With this choice of bases, sparsity in the representation of many data types is 
to be generally expected.^ 

In Figure [3] we explore coefficient statistics and the effects of thresholding on the reconstruction errors for 
a model of terrestrial topography. We focus on the fifth, or "North- American" chunk of our cubed sphere, and 
use the orthogonal D2 (Haar), D4 and D6 wavelet base^^on the interval, with preconditioning.^^ The top row 
uses the common conventions in plotting the wavelet and scaling coefficients in each of the bases after (hard) 
thresholding^ such that only the coefficients larger than their value at the 85*^ percentile level survive. The 
coefficients that have now effectively been zeroed out are left white in these top three panels. The middle series 
of panels of Figure [3] plots the spatial reconstruction after thresholding at this level; the root mean squared (rms) 
errors of these reconstructions are quoted as a percentage of the original root mean squared signal strengths. 
The thresholded wavelet transforms allow us to discard, as in these examples, 85% of the numbers required to 
make a map of North American topography in the cubed-sphere pixel basis: the percentage error committed is 
only 6.3%, 5.2% and 9.2% according to this energy criterion in the D2, D4 and D6 bases, respectively. From 
the map views it is clear that despite the relatively small error, the D2 basis leads to unsightly block artifacts in 
the reconstruction, which are largely avoided in the smoother and more oscillatory D4 and D6 bases. A view of 
the coefficient statistics is presented in the lowermost three panels of Figure [3] The coefficients are roughly log- 
normally distributed, which helps explain the success of the thresholded reconstruction approach. We conclude 
that the D4 basis is a good candidate for geophysical data representation. 

For applications in seismology we now illustrate the performance of the cubed-sphere wavelet basis in com- 
pressing seismological Earth models such as the one by Montelli et al. (2006) of compressional (P) wavespeed 
heterogeneity and another by Ritsema et al. (2010) of shear (5) wavespeed perturbations. At a depth of about 
670 km. Figs [4^ and|4^ show the wavespeed anomalies from the average at that depth. The wavelet transform in 
the D4 basis (with special boundary filters and after preconditioning, and up until scale J = 3) was thresholded 
and the results re-expanded to the spatial grid, identically as we did for the topography in Figure [3) The results 
for specific values of the thresholding (quoted as the percentile of the original wavelet coefficients) are shown 



in Figs|4]3 andjZjF for the 85*^, and FigsjZ]^ and|4^ for the 95*^, respectively. At each level of thresholding the 
number of nonzero wavelet /scaling expansion coefficients is quoted: at 0% thresholding this number is identical 
to the number of pixels in the surficial cubed sphere being plotted. 

The reconstruction error can be visually assessed from the pictures; it is also quoted next to each panel as 
the percentage of the root mean squared error between the original and the reconstruction, normalized by the 
root mean squared value of the original in the original pixel representation, in percent. Specifically, we calculate 
and quote the ratio of ^2 norms in the pixel-basis model vector m, 

100 X ||m-5{r[A(m)]}||2/||m||2, (2) 

which, in the lower-right annotations is called the "% error norm" . We write A for any of the wavelet (analysis) 
transforms that can be used and S (synthesis) for their inverses, and T for the hard thresholding of the wavelet 
and scaling coefficients. In Figs|4^ and[^, the misfit quantity ([2| is represented as a black line relevant to the 
left ordinate labeled "^2 error norm", which shows its behavior at 1% intervals of thresholding; the filled black 
circles correspond to the special cases shown in the map view. Only after about 80% of the coefficients have 
been thresholded does the error rise above single-digit percentage levels, but after that, the degradation is swift 
and inexorable. The blue curves in Figs|4^ and|^ show another measure relevant in this context, namely the 
ratio of the £1 norms of the thresholded wavelet coefficients compared to the original ones, in percent, or 

100 X \\T[A{ui)]\y\\A{xn)\\,. (3) 

The £2 ratios ([2| in the black curves (and the left ordinate) evolve roughly symmetrically to the ii ratios ([3| 
in the blue curves (and the right ordinate), though evidently their range is different. The third measure that is 
being plotted as the red curve is the "total variation" norm ratio, in percent, namely 

100 X ||V5{r[^(m)]}||,/||Vm||,, (4) 

whereby ||Vm||-|^ is the sum over all voxels of the length of the local gradient of m. By this measure, which 
is popular in image restoration applications,'^^'^^ the quality of the reconstruction stays very high even at very 
elevated levels of thresholding; we note that its behavior is not monotonic and may exceed 100%. 

Finally, the new wavelet construction can be used to study the joint properties of the wavelet coefficients 
of seismic wavespeed models, as we illustrate in Figure [5] There, we report the correlation between wavelet 
coefficients in the Montelli and Ritsema models as a function of scale and approximate geographical position 
(see again Figure |2] for the numbering scheme of the cubed-sphere chunk). A rendering of the two-dimensional 
density of the data is accompanied by the value of their correlation coefficient (lower left labels) where this is 
deemed significant at the 95% level, and the slope of the total-least-squares based fit in this space (upper right 
labels), which is only quoted when the correlation coefficients exceeded 0.35. This should provide an estimate 
of the logarithmic ratio of shear- wave to compressional-wave speed perturbations, dlnVs/SlnVp^ an important 
discriminant in the interpretation of the (thermal or chemical) cause of seismic velocity anomalies. The variation 
of this ratio as a function of scale and chunk position yields information that will be of use for geochemical and 
geodynamical studies, and the orthogonality of the wavelet basis in scale and physical space removes some 
of the arbitrariness in the calculation. Ritsema's model does not have much structure at the smallest scales. 
From scale 3 onward a positively correlated pattern begins to emerge, though even at this particular scale, the 
correlation coefficients remain below the relatively stringent 0.35 level. Wavelets and scaling coefficients are well 
correlated at the largest scale 4 considered, with several of the correlation coefficients exceeding our threshold. 
This information is useful to geophysicists,^^ especially given our new-found ability to study the regional 
variation of such ratios, taking into account their dependence on scale length. 

Simons et al. (2011) detail the reasoning behind the construction of the "superchunk" cubed sphere shown 
in Figure [l] (n^/it). In a nutshell, the cubed-sphere bases are to be used not simply for the representation and 
analysis of seismic models, but also to parameterize the inversion of primary data for such models. Since the 
edge-cognizant transforms are not norm-preserving, the thresholding steps involved in using common algorithms 
such as (F)ISTA)llI3ll2l if unmodified, lead to artifacts in the solution which are largely avoided by building the 
wavelet transforms in the superchunk domain and ignoring the edges altogether.— 



3. THE SPHERICAL SLEPIAN TREE TRANSFORM 



Geophysical and cosmological signals are constrained by the physical processes that generate them and con- 
ditioned by the sensors that observe them. On the real line and in the plane, both physical and sampling 
constraints lead to assumptions of a bandlimit: that a signal contains no energy outside some supporting region 
in the frequency domain. Bandlimited signals on the sphere are zero save for the low-frequency spherical- 
harmonic components. Spherical harmonics are not orthogonal on arbitrary portions of the sphere, yet for the 
study of geophysical processes we may not have access to nor interest in signal originating from outside a specific 
geographic region of interest. By construction, Slepian functions satisfy the dual constraints of bandlimitation 
and spatial concentration by optimization of an energy criterion that concentrates as much energy as possible 
into the region of interest for a given bandlimit .^^^ With Slepian function bases, however, it is "all or nothing": 
given a certain bandwidth and a spatial region of interest, the functions (variably) fill the entire bandwidth range 
and ultimately the entire spatial target. The wavelet constructions that we introduced in Section |2j on the other 
hand, were compactly supported and had multiresolution properties which the Slepian basis does not possess. 

Here we develop an algorithm for the construction of Slepianesque dictionary elements that are bandlimited, 
localized, and multiscale. It is based on a subdivision scheme that constructs a binary tree from subdivisions of 
the region of interest. Such a dictionary has many nice properties: it closely overlaps with the most concentrated 
Slepian functions on the region of interest, and most element pairs have low coherence. Therefore, they, too, 
should be eminently suitable to solve ill-posed inverse problems in geophysics and cosmology. Though the new 
dictionary is no longer composed of purely orthogonal elements like the Slepian basis, it can also be combined 
with modern inversion techniques that promote sparsity in the solution. Thereby they provide significantly lower 
residual error after reconstruction when compared to "classically optimal" Slepian inversion techniques.*^ 

3.1. A Multiscale Dictionary of Slepian Functions 

We now turn our focus to numerically constructing a dictionary V of functions that can be used to approximate 
bandlimited signals on the sphere and allows for the reconstruction of signals from their point samples. Let 
7^ C 5^ be a simply connected subset of the sphere. With L the bandwidth, the dictionary V will be composed 
of functions bandlimited to spherical harmonic degrees < / < L. The construction is based on a binary tree. 
The node capacity, n, is a positive integer. Each node of the tree corresponds to the first n Slepian functions 
with bandlimit L and concentrated on a subset IZ' C IZ. The top tree node is for the entire region 7^, and each 
node's children correspond to a division of IZ' into two roughly equally sized subregions. As the child nodes 
will be concentrated in disjoint subsets of IZ' ^ all of their corresponding functions and children are effectively 
incoherent. We now fix a height H of the tree: the number of times to subdivide IZ. The height is determined 
as the maximum number of binary subdivisions of IZ that can have n well concentrated functions. That is, we 
find the minimum integer H such that n > N2-h\^\^i^., with N the Shannon number, which has the solution 
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A complete binary tree with height H has 2^+^ — 1 nodes, so from now on we will denote the dictionary 

as the set of |2^7^, L,n I = ^ (2^+^ — 1) functions thus constructed on region IZ with bandlimit L and node capacity 
n. Figure |6] shows the tree diagram of the subdivision scheme. We use the standard enumeration of nodes wherein 
node (j, •) is subdivided into child nodes (2j, •) and (2j + 1, •), and at a level < h < the nodes are indexed 
from 2^ <j < 2^+^ - 1. More specifically, for j = 1, 2, . . ., we have U^^^ = VS'^^^ U 7^(2^■+^). Furthermore, letting 
be the a^^ Slepian function on IZ' (the solution to the classical Slepian concentration criterion eq. 4.1 in 



ref. 13, with concentration region IZ')^ we have that 

Figure [t] shows an example of the construction when IZ is the African continent. Note how, for example, d^^'^^ 
and d!^^'^^ are the first Slepian functions associated with the subdivided domains of VS^\ 



To complete the top-down construction, it remains to decide how to subdivide a region VJ into equally sized 
subregions. For roughly circular connected domains, the first Slepian function has no sign changes, and the 
second Slepian function has a single zero-level curve that subdivides the region into approximately equal areas; 
when VJ is a spherical cap, the subdivision is exact We thus subdivide a region VJ into the two nodal domains 
associated with the second Slepian function on that domain; see Figure [S] for a visualization of the subdivision 
scheme as applied to the African continent. 

3.2. Concentration, Range, and Incoherence 

The utility of the tree construction presented above depends on its ability to represent bandlimited functions 
in a region 7^, and its efficacy at reconstructing functions from point samples in IZ. These properties, in turn, 
reduce to questions of concentration, range, and incoherence. First, dictionary V is concentrated in IZ if its 
functions are concentrated in IZ. Second, the range of dictionary V is the subspace spanned by its elements. 
Ideally, the basis formed by the first few Slepian functions on 7^ is a subspace of the range of V. Third, when V 
is incoherent, pairwise inner products of its elements have low amplitude: pairs of functions are approximately 
orthogonal. This, in turn, is a useful property when using V to estimate signals from point samples. 

Unlike the concentration eigenvalues corresponding to the Slepian functions on 7^, not all of the eigenvalues of 
the elements of reflect their concentration within this top-level (parent) region. We thus define the modified 
concentration value 

zy^^>) = / U^^^'^Xx)]^ d^{x). (5) 



Recalling that II2 = "^^i^ value v is simply the percentage of energy of the (j, a)^^ element that is 

concentrated in IZ. This value is always larger than the element's eigenvalue, which relates its fractional energy 
within the smaller subset VS^\ 

The size of dictionary Vn,L,n is generally larger than the Shannon number A^|7^|,L = {L ^ 1)^/ \1Z\ /4/7r for 
any node capacity n, and as a result it cannot form a proper basis: it has too many functions. Ideally, then, we 
require that elements of the dictionary span the space of the first N\^\^l Slepian functions. One possible answer 
to the question of the range of V is given by studying the angle between the subspaces^^ spanned by elements 
of V and the first a functions of the Slepian basis, for a = 1, 2, — The angle between two subspaces A and 5, 
possibly with different dimensions, is given by the formula 

/(A, 5) = min Tsup Z(x, 5), sup Z(?/, A)^ , whereZ(x,5)= inf Z(x, y) = cos""*- n^^^^ • (6) 
\xeA yeB ) y^B \\x\\ 

Here, Pb is the orthogonal projection onto space B and all of the norms are with respect to the given subspace. 
The angle Z(A,5) is symmetric, non- negative, and zero \S. A (Z B or B (Z A] furthermore it is invariant under 
unitary transforms applied to both on A and 5, and admits a triangle inequality. It is thus a good indicator of 
distance between two subspaces; furthermore, it can be calculated accurately given two matrices whose columns 
span A and B. We can therefore identify the matrices A and B with the subspaces spanned by their columns. 

Let {Gn,L)i.^ denote the matrix containing the first a column vectors which are the spherical harmonic 
expansion coefficients of the traditional Slepian functions for a given region IZ and a bandwidth L. Let D denote 
the (L + l)^x |P7^, L,n[ matrix containing the spherical harmonic representations of the elements of the dictio- 
nary P7^,L,n• Figure 9 shows l{Gi:a^D) for 7Z = Africa with L = 36. The Shannon number is ^Africa 36 ^ 
It is clear that while the dictionaries T^Africa 36 1 -^Africa 36 2 strictly span the space of functions 

bandlimited to L = 36 and optimally concentrated in Africa, they are a close approximation: the column span 
(^Africa, 36) l a nearly linearly dependent with the spans of ^Africa,36,i ^Africa, 36, 2' ^ significantly 
larger than' 7VAfrica,36- 

The requirement that the dictionary elements form an orthogonal basis is less important than the property of 
mutual incoherence. This we identify when the inner product between pairs of elements is almost always very low. 



Figure 10 shows that the two tree constructions on continental Africa have good incoherency properties: most 
dictionary element pairs are nearly orthogonal. As can be seen from Figure [TO^, most pairwise inner products are 
nearly zero. More specifically, as expected, dictionary elements (j, 1), (2j, 1), (2j + 1, 1), (2(2j), 1), (2(2j) + 1, 1), 



(2(2j + 1), 1), (2(2j + 1) + 1, 1), . . ., tend to have large inner products, while those with non-overlapping borders 
do not. This exact property is also visible in the two diagonal submatrices of Figure [ToJd. In the off-diagonals, 
due to the orthonormality of the construction, elements of the form (j, 1) and (j, 2) are orthogonal. In contrast, 
due to the nature of the tree subdivision scheme, elements of the form (2j, 1) or (2j + 1,1) and (j, 2) have a large 
inner product. However, the number of connections between nodes and their ancestors is 0{n [2^i7]), while the 
total number of pairwise inner products is (9([n2^]^); and for reasonably sized values of L the ratio of ancestral 
connections to pairwise inner products gets to be small, see Figure [TT] 

4. CONCLUSIONS 

We have introduced two classes of spherical parameterizations and discussed their properties in the context of 
geophysical model analysis and representation. The first involved the porting of traditional Cartesian wavelet 
transforms to the sphere via a cubed-sphere mapping, with and without special consideration for the boundaries 
between the six chunks constituting the entire spherical surface. The second was a novel elaboration of the 
classical ideas of signal concentration on the sphere. Starting from optimally spatially concentrated bandlimited 
spherical "Slepian" functions, we construct a dictionary of functions occupying a binary tree, where the elements 
of the tree are successive levels of Slepian functions calculated for an increasingly subdivided spatial domain. 

Both the wavelet bases and the Slepian tree frames that we discuss in this paper are flexible and efficient 
ways of studying the information content and spatiospectral structure of geophysical and cosmological data. But 
in both cases their utility will also be derived from using them in the parameterization of geophysical inverse 
problems. In ref. [l]we discuss an iterative algorithm that solves an inverse problem in seismic tomography while 
promoting sparsity in the cubed-sphere wavelet basis in which the unknown model is expressed. Ref . |2] explores 
the inverse problem of approximating bandlimited, heterogeneously concentrated, essentially multiscale functions 
on the sphere from their samples. 

In this paper and elsewhere, we have shown by example that many geophysical signals are sparse in the wavelet 
and Slepian (tree) "bases". To tackle large-scale and potentially ill-conditioned inverse problems, the sparsity 
of the solution should actively be encouraged. In either case the inverse problem comes down to minimizing a 
mixed £2-^1 functional of the form 

^(w) = ||K-S-w-d||i + 2A||w||i, (7) 

whereby d is a data set, K a linear operator, and S a certain synthesis map that transforms the set of unknown 
coefficients w into the domain of K. Thus, S could relate to the wavelet bases introduced in Section [2] or to the 
Slepian-tree dictionary introduced in Section [3j Numerical experiments conducted for a variety of experimental 
setups under geophysically realistic conditions have already indicated that both types of constructions will be 
amenable to solving inverse problems in geophysics and beyond. Further research will be reported elsewhere. 
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Figure 1. Aerial view showing our adaptations of the cubed sphere of ref . [5] Of the front-facing four of the in total six 
"chunks", one is gridded to reveal its 2^^ surface elements (N = 4). The thick black lines identify the boundaries of the 
six "chunks". The thin gray lines correspond to the boundaries of the overlapping "superchunks" as discussed in the text. 




Figure 2. Geometry, nomenclature. 



and 



numbering of the six faces of our cubed sphere in a two-dimensional view. 
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Figure 3. Wavelet and scaling coefficients (top), space- domain reconstructions after thresholding (middle), and "signed" 
histograms (bottom) of the wavelet and scaling coefficients of the "North American" face of the cubed-sphere version of 
the Earth's topography, to a 2'^ dyadic subdivision with J = 3. We used preconditioned interval wavelet transforms. 
All coefficients were hard-thresholded, retaining only the 15% largest coefficients by absolute value. In the top row, the 
locations of zeroed coefficients are rendered white; those are also captured by the white bars in the histograms. The root 
mean squared (rms) error of the reconstruction after thresholding is indicated as a percentage of the signal rms. Tick 
marks on the color bars identify the 5*^, 50^^ and 95*^ percentile of the coefficients or the spatial reconstructions after 
thresholding, respectively. Interior ticks on the histograms roughly coincide with these same percentiles as applied to 
either the positive and negative coefficients when expressed on a logarithmic scale. Histograms for the positive coefficients 
point up and have ordinates in positive percentages, histograms for the negative coefficients point down and have ordinates 
in negative percentages; these percentages are with respect to the total number of positive and negative coefficients. The 
blue and red shaded areas of the histograms reffect the coefficients retained at the global 85*^ thresholding level. 
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Figure 4. Sparsity and reconstruction stability of two global seismic wavespeed models under incre menta l hard thresh- 
olding of their wavelet and scaling coefficients using the preconditioned edge-cognizant D4 wavelet basi£3^in the angular 
coordinates of the cubed sphere, as developed in this paper, (a-d) Results for the compressional wave speed seismic model 
of ref . [43I and (e-h) for the shear wave speed seismic model of ref. at the same depth of 677 km below the surface of 
the Earth, for cubed spheres with 6 x 2^^ elements (A^ = 7), and to a 2'^ dyadic subdivision (J = 3). As a function of the 
percentage of the coefficients that are being thresholded, and relatively to the original unthresholded values, the bottom 
panels quote the spatial £2 norms of the reconstruction error (in black), the total variation norms of the reconstructed 
images in the space domain (in red), and the ii norms of the coefficients that remain (in blue). The reconstructions 
remain faithful to the originals even at elevated levels of thresholding. 
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Figure 5. Joint properties of seismic mantle structure in the Montelli (2006) P-wave and Ritsema (2010) S'-wave speed 
models, at 677 km depth in the Earth. Every row corresponds to a different scale in the D4 wavelet decomposition of the 
models. Each panel shows the logarithmic density of observations. Black shading corresponds to the maximum density 
in each panel; all patches that account for less than 1% of the observations are rendered white. Using total least squares 
a regression line was fit to all sets with correlation coefficients exceeding 0.35. The slope of the line, a measure of the 
SlnVs /SlnVp ratio appears in the top right corners. Significant correlation coefficients are quoted in the bottom left. 
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Figure 6. The binary tree subdivision scheme and associated dictionary T>tz,l for the multiscale Slepian tree transform. 
We define the top-level region IZ as IZ^-^^ and the generic subsets IZ' as 
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Figure 7. Slepian tree dictionary ^Africa 36 i ^^^^ 1^1 ~ showing the functions d*^^'^-* through d*^^'^-* and d'^^^^'^-* 
through The X-axis is longitude, the y-axis is colatitude. Regions of concentration jT^*^*-*! are outhned. 
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Figure 8. Second {a — 2) Slepian functions associated with the regions in Figure [7^-f. Regions of concentration, and the 
dividing contour (the zero-level set), are drawn in green. Blue and red represent the sign of the Slepian function values. 
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Figure 9. Angles (in degrees) between the spaces spanned by {G\'fi\^3^ ^ and the dictionary matrices D-jz^se^i {left) 
and -D7j,36,2 (right), TZ = Africa, a = 1, . . . jSN^-pii^L- Thick Unas correspond to integer multiples of the Shannon num- 
ber ^Africa,36 ~ 79. 




Figure 10. Magnitudes of pairwise inner products of dictionaries (a) ^Africa 36 i -^Africa 36 2' (^)' ^^i^k 

lines separate inner products between the 127 elements with a = 1 (top left) and a = 2 (bottom right), and their cross 
products. Values are between (white) and 1 (black). 
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Figure 11. Empirical cumulative distribution functions of pairwise inner product magnitudes for dictionaries ^Africa 36 i 
(left) and ^Africa 36 2 i^'^d^t). In both cases, approximately 95% of pairwise inner products have a value within ±0.1. 



